home *** CD-ROM | disk | FTP | other *** search
/ CD ROM Paradise Collection 4 / CD ROM Paradise Collection 4 1995 Nov.iso / science / fastmap.zip / FASTMAP.C < prev    next >
Text File  |  1992-09-14  |  36KB  |  1,186 lines
  1. /* FASTMAP.C */
  2. /* Copyright (C) David Curtis 1992 */
  3.    
  4. #ifndef FASTMAP_H
  5. #include "fastmap.h"
  6. #endif
  7.  
  8. #include <math.h>
  9.  
  10. #include <string.h>
  11.  
  12. #if 0
  13. #define strchr index
  14. /* use this if your library contains index() and not strchr() */
  15. #endif
  16.  
  17. #if __ZORTECH__
  18. unsigned _stack=30000;
  19. #endif
  20.  
  21. #if __ZORTECH__  || __STDC__
  22. #include <stdlib.h>
  23. #else
  24. extern char *calloc PROTO((unsigned nelem,unsigned elsize));
  25. #endif
  26.  
  27. #ifndef max
  28. #define max(x,y) ((x)>(y)?(x):(y))
  29. #define min(x,y) ((x)<(y)?(x):(y))
  30. #endif
  31.  
  32. #define BIG 1e8
  33.  
  34. #if USE_ENUMS
  35. enum result_index { rr,nn,nr,rn,ru,nu,ur,un,uu };
  36. #else
  37. #define rr 0
  38. #define nn 1
  39. #define nr 2
  40. #define rn 3
  41. #define ru 4
  42. #define nu 5
  43. #define ur 6
  44. #define un 7
  45. #define uu 8
  46. #endif
  47.  
  48. float fixdist[MAXDISTS];
  49. int nfixdist;
  50. struct marker_t disease,*marker;
  51. FILE *input_file,*output_file,*graph_file,*debug_file;
  52. int num_points,num_peds;
  53. float mind,maxd;
  54. float **tablods,*totlods;
  55. float total_meioses[MAXPEDS];
  56. float reliability[MAXPEDS];
  57.  
  58. #define GRATIO 0.61803399
  59. #define MGRATIO (1.0-GRATIO)
  60. #define SHFT(a,b,c,d) (a)=(b);(b)=(c);(c)=(d);
  61. /* GRATIO, MRATIO and SHFT are used for the golden ratio minimisations.
  62. These are based on the routine described in Numerical Recipes in C by
  63. Press, Flannery, Teukolsky and Vettering, Cambridge University Press. 
  64. */
  65.  
  66. #define error(s,t) error_func(__LINE__,__FILE__,s,t)
  67.  
  68. void print_consts()
  69. {
  70. printf("FASTMAP - program constants:\n\n\
  71. MAXPEDS      = %d\n\
  72. MAXMARKERS   = %d\n\
  73. MAXPAIRS     = %d\n\
  74. MAXDISTS     = %d\n\
  75. MINFRACTION  = %6.4f\n\n\n",
  76. MAXPEDS, MAXMARKERS, MAXPAIRS, MAXDISTS, MINFRACTION);
  77. }
  78.  
  79. /* allocate memory for marker data */
  80.  
  81. int init_marker(m,numf)
  82. struct marker_t *m;
  83. int numf;
  84. {
  85. int i;
  86. m->rec=m->nonrec=0;
  87. m->lods=0;
  88. if ((m->npairs=(int *)calloc(numf,sizeof(int)))==0) return 0;
  89. if ((m->rec=(float *)calloc(numf,sizeof(float)))==0) return 0;
  90. if ((m->nonrec=(float *)calloc(numf,sizeof(float)))==0)
  91.     { free(m->rec); m->rec=0; return 0; }
  92. if ((m->lods=(struct lod_rf_t **)calloc(numf,sizeof(struct lod_rf_t *)))==0)
  93.     { free(m->nonrec); free(m->rec); m->nonrec=m->rec=0; return 0; }
  94. for (i=0;i<numf;++i)
  95.    if ((m->lods[i]=(struct lod_rf_t *)calloc(MAXPAIRS,sizeof(struct lod_rf_t)))==0)
  96.      {
  97.      while (--i>=0) free(m->lods[i]);
  98.      free(m->lods);
  99.      m->lods=0;
  100.      free(m->nonrec);
  101.      free(m->rec);
  102.      m->nonrec=m->rec=0;
  103.      return 0;
  104.      }
  105. return 1;
  106. }
  107.  
  108.  
  109. void clear_2p(m)
  110. info2p_t *m;
  111. {
  112. m->rec=m->nonrec=m->fraction_used=0.0;
  113. }
  114.  
  115.  
  116. /* next two functions set up likelihoods for two-points and three-
  117. points with recombination fraction of 50% with disease locus, to avoid 
  118. recalculating for each lod score required */
  119.  
  120. void set_null_like_2p(m)
  121. info2p_t *m;
  122. {
  123. if (!m->rec && !m->nonrec) m->null_like=0.0;
  124. else m->null_like=log_bin_like(m->rec,m->rec+m->nonrec,0.5);
  125. }
  126.  
  127. void set_null_like_3p(mm)
  128. info3p_t *mm;
  129. {
  130. float theta;
  131. if (!mm->rrnn[rr] && !mm->rrnn[nn] && !mm->rrnn[nr] && !mm->rrnn[rn])
  132.   mm->null_like=0.0;
  133. else
  134.   {
  135.   theta=anti_kosambi(fabs(mm->l_position-mm->r_position));
  136.   mm->null_like=(mm->rrnn[rr]+mm->rrnn[nn])*log10(0.5*(1-theta))+
  137.                 (mm->rrnn[nr]+mm->rrnn[rn])*log10(0.5*theta);
  138.   }
  139. }
  140.  
  141. /* calculate three-point lod score with disease at recombination phi 
  142. with left marker and pi with right marker */
  143.  
  144. float lod_3p(d,mm,rel)
  145. struct marker_t *d;
  146. info3p_t *mm;
  147. double rel;
  148. {
  149. float lod,phi,pi;
  150. if (!mm->rrnn[rr] && !mm->rrnn[nn] && !mm->rrnn[nr] && !mm->rrnn[rn])
  151.   return 0.0;
  152. phi=anti_kosambi(fabs(d->position-mm->l_position));
  153. pi=anti_kosambi(fabs(d->position-mm->r_position));
  154. lod=mm->rrnn[rr]*log10(phi*pi*rel+(1.0-pi)*(1.0-phi)*(1.0-rel));
  155. lod+=mm->rrnn[nr]*log10((1-phi)*pi*rel+(1.0-pi)*phi*(1.0-rel));
  156. lod+=mm->rrnn[rn]*log10(phi*(1-pi)*rel+pi*(1-phi)*(1.0-rel));
  157. lod+=mm->rrnn[nn]*log10((1-phi)*(1-pi)*rel+pi*phi*(1.0-rel));
  158.  
  159. lod-=mm->null_like;
  160. return lod;
  161. }
  162.  
  163. /* calculate two-point lod score with disease at recombination theta 
  164. with marker */ 
  165.  
  166. float lod_2p(d,m,rel)
  167. struct marker_t *d;
  168. info2p_t *m;
  169. double rel;
  170. {
  171. float lod,theta;
  172. if (!m->rec && !m->nonrec) return 0.0;
  173. theta=anti_kosambi(fabs(d->position-m->position));
  174. lod=log_bin_like(m->rec,m->rec+m->nonrec,theta*rel+(1.0-theta)*(1.0-rel));
  175. lod-=m->null_like;
  176. return lod;
  177. }
  178.  
  179. /* recombination fraction between the two markers referred to by info 
  180. structures r and l */
  181.  
  182. float inf_rf_between(r,l)
  183. info2p_t *r,*l;
  184. {
  185. return anti_kosambi(fabs(r->position-l->position));
  186. }
  187.  
  188. /* calculate PIC value from allele frequencies */
  189.  
  190. float get_prob_inf(p,n)
  191. float p[];
  192. int n;
  193. {
  194. int i,j;
  195. float pic;
  196. pic=1;
  197. for (i=0;i<n;++i)
  198.   {
  199.   pic-=p[i]*p[i];
  200.   }
  201. for (i=0;i<n-1;++i)
  202.   for (j=i+1;j<n;++j)
  203.     {
  204.     pic-=2*p[i]*p[i]*p[j]*p[j];
  205.     }
  206. return pic;
  207. }
  208.  
  209. /* get recombination fraction corresponding to distance in centimorgans
  210. */
  211.  
  212. float anti_kosambi(dist)
  213. double dist;
  214. {
  215. float ed;
  216. ed=exp(dist/25.0);
  217. return 0.5*(ed-1.0)/(ed+1.0);
  218. }
  219.  
  220. /* routine to get one line of input */
  221.  
  222. char *get_line(line,maxlen)
  223. char *line;
  224. int maxlen;
  225. {
  226. char buff[1001];
  227. *buff='\0';
  228. fgets(buff,1001,input_file);
  229. strncpy(line,buff,maxlen);
  230. line[maxlen]='\0';
  231. if (strchr(line,'\n')) *strchr(line,'\n')='\0';
  232. return line;
  233. }
  234.  
  235. input_marker(m,num_peds)
  236. struct marker_t *m;
  237. int num_peds;
  238. {
  239. char buff[241],line[241],rest[241];
  240. float freq[100];
  241. int nargs,nall,i;
  242. prompt("\
  243. Enter marker name, position and heterogeneity value or allele frequencies\n\
  244. (blank line all markers entered)\n");
  245. get_line(line,240);
  246. nargs=sscanf(line,"%20s %f %[^\n]",m->name,&m->position,rest);
  247. if (nargs<1) return 0;
  248. else if (nargs<3)
  249.    return error("Bad argument format:\n",line);
  250. else if (!init_marker(m,num_peds))
  251.    return error("Not enough memory for this marker:\n",line);
  252. else 
  253.   {
  254.   float tot;
  255.   nall=0;
  256.   tot=0.0;
  257.   do {
  258.   strcpy(buff,rest);
  259.   if ((nargs=sscanf(buff,"%f %[^\n]",&freq[nall],rest))<1) break;
  260.   if (freq[nall]>1.0)
  261.      return error("Heterogeneity or allele frequency greater than 1:\n",line);
  262.   tot+=freq[nall++];
  263.   }  while (nargs==2);
  264.   if (nall==0) return error("No heterogeneity value:\n",line);
  265.   else if (nall>1 && tot!=1.0) return error("Allele frequencies do not sum to 1:\n",line);
  266.   else if (nall==1)  m->prob_inf=freq[0];
  267.   else m->prob_inf=get_prob_inf(freq,nall);
  268.  
  269. if (debug_file) fprintf(debug_file,"%s.prob_inf=%f\n",m->name,m->prob_inf);  
  270.  
  271.   for (i=0;i<num_peds;++i)
  272.     {
  273.     float lod,rf,maxlod;
  274.     prompt("Enter a number of theta/lod score pairs or theta and max lod:\n\
  275. (blank line if genotyping not done in this family):\n");
  276.     get_line(line,240);
  277.     m->npairs[i]=0;
  278.     if (sscanf(line," %f %f",&rf,&lod)<1)
  279.        continue; /* genotyping not done in this family */
  280.     maxlod=0;
  281.     strcpy(rest,line);
  282.     do {
  283.     strcpy(buff,rest);
  284.     *rest='\0';
  285.     if ((nargs=sscanf(buff,"%f %f %[^\n]",&rf,&lod,rest))<2) break;
  286.     m->lods[i][m->npairs[i]].rf=rf;
  287.     m->lods[i][m->npairs[i]].lod=lod;
  288.     ++m->npairs[i];
  289.     maxlod=max(fabs(lod),maxlod);
  290.     } while (nargs>=2);
  291.     if (maxlod<0.01)  /* marker uninformative */
  292.       {
  293.       m->lods[i][0].rf=0;
  294.       m->lods[i][0].lod=0;
  295.       m->npairs[i]=1;
  296.       }
  297.     }
  298.   }
  299. return 1;
  300. }
  301.  
  302. /* estimate the equivalent number of nonrecombinant and recombinant 
  303. phase-known meoises which would approximate to the input lod scores, 
  304. given a supplied value for "reliability" with which phenotype 
  305. corresponds with genotype */
  306.  
  307. float get_nr_from_lods(m,ped,rel)
  308. struct marker_t *m;
  309. int ped;
  310. double rel;
  311. {
  312. float N,elod,ssq;
  313. ssq=0.0;
  314. if (m->npairs[ped]==0) ; /* marker not tested in this pedigree */
  315. else if (m->npairs[ped]==1)  /* only max lod/theta */
  316.   {
  317.   if (m->lods[ped][0].lod==0.0 && m->lods[ped][0].rf==0.0)
  318.     m->rec[ped]=m->nonrec[ped]=0.0; 
  319.     /* marker uninformative in this pedigree */
  320.   else if (m->lods[ped][0].lod<=0.0)
  321.     error("If only one lod it must be positive","");
  322.   else
  323.     {
  324.     elod=log_bin_like(m->lods[ped][0].rf*rel+
  325.              (1.0-m->lods[ped][0].rf)*(1.0-rel),
  326.              1,(m->lods[ped][0].rf*rel+(1-m->lods[ped][0].rf)*(1-rel)));
  327.     elod-=log_bin_like(m->lods[ped][0].rf*rel+
  328.              (1.0-m->lods[ped][0].rf)*(1.0-rel),1,0.5);
  329. /* calculate what lod score would be with only one informative 
  330. meiosis */
  331.     N=m->lods[ped][0].lod/elod;
  332.     m->rec[ped]=N*(m->lods[ped][0].rf*rel+
  333.                      (1.0-m->lods[ped][0].rf)*(1.0-rel));
  334.     m->nonrec[ped]=N*((1.0-m->lods[ped][0].rf)*rel+
  335.                      m->lods[ped][0].rf*(1.0-rel));
  336. /* then get total number of informative meioses and divide into 
  337. recombinants and nonrecombinants */
  338.     }
  339.   }
  340. else if (m->npairs[ped]==2)
  341.   {
  342.   float a,b,c,d;
  343.   a=log10(m->lods[ped][0].rf*rel+(1-m->lods[ped][0].rf)*(1-rel))-log10(0.5);
  344.   c=log10(m->lods[ped][1].rf*rel+(1-m->lods[ped][1].rf)*(1-rel))-log10(0.5);
  345.   b=log10((1-m->lods[ped][0].rf)*rel+m->lods[ped][0].rf*(1-rel))-log10(0.5);
  346.   d=log10((1-m->lods[ped][1].rf)*rel+m->lods[ped][1].rf*(1-rel))-log10(0.5);
  347.   m->rec[ped]=(d*m->lods[ped][0].lod-b*m->lods[ped][1].lod)/(a*d-b*c);
  348.   m->nonrec[ped]=(a*m->lods[ped][1].lod-c*m->lods[ped][0].lod)/(a*d-b*c);
  349. /* if two pairs of values supplied then solve the simultaneous equation
  350. to produce an exact solution for the numbers of meioses which would 
  351. give these values */
  352. /* it is probably best to avoid using this routine, as the curve seems 
  353. likely to be wildly inaccurate at other recombination fractions */
  354.   }
  355. else
  356.   {
  357. /* use a numerical method to obtain the likely proportion of 
  358. recombinants, then calculate the total number of meioses from that */
  359.   float ratio;
  360.   ratio=get_rf_fraction(m->lods[ped],m->npairs[ped],rel,&ssq,0.01);
  361.   N=get_num_meioses(m->lods[ped],m->npairs[ped],rel,ratio);
  362.   ssq*=N*N;
  363.   m->rec[ped]=N*ratio;
  364.   m->nonrec[ped]=N*(1-ratio);
  365.   }
  366. return ssq;
  367. /* if the last method is used, return a measure of the closeness of fit
  368. obtained to the input values */
  369. }
  370.  
  371. /* use the above routine to select the best fitting reliability value 
  372. (over all markers) unless it is already supplied, and using this to
  373. settle on the overall best estimates for the numbers of nonrecombinant
  374. and recombinant meioses for each marker */ 
  375.  
  376. meioses_from_lods(numf,num_markers)
  377. int numf,num_markers;
  378. {
  379. int i,j;
  380. for (i=0;i<numf;++i)
  381.   {
  382.   if (reliability[i]>=0)
  383. /* reliability value supplied for this pedigree */
  384.     for (j=0;j<num_markers;++j)
  385.        get_nr_from_lods(&marker[j],i,reliability[i]);
  386.   else
  387.     {
  388. /* need to fit value across all markers */
  389.     float f,f0,f1,f2,f3,x0,x1,x2,x3,fmin;
  390.     float minrel;
  391.     int canfit;
  392.     minrel=1.0;
  393.     for (canfit=0,j=0;j<num_markers;++j)
  394.       if (marker[j].npairs[i]>2) { canfit=1; break; }
  395.     if (canfit) /* at least one marker with >=3 pairs */
  396.       {
  397.       /* use a golden ratio minimisation */
  398.       x0=0.5;
  399.       x3=1.0;
  400.       x1=0.5+0.5*MGRATIO;
  401.       x2=0.5+0.5*GRATIO;
  402.       for (f1=0.0,j=0;j<num_markers;++j)
  403.           if (marker[j].npairs[i]>2) f1+=get_nr_from_lods(&marker[j],i,x1);
  404.       for (f2=0.0,j=0;j<num_markers;++j)
  405.           if (marker[j].npairs[i]>2) f2+=get_nr_from_lods(&marker[j],i,x2);
  406.       while (0.005<fabs(x1-x2))
  407.         {
  408.         if (f2 < f1)
  409.           {
  410.           SHFT(x0,x1,x2,GRATIO*x1+MGRATIO*x3)
  411.           for (f=0.0,j=0;j<num_markers;++j)
  412.              if (marker[j].npairs[i]>2) f+=get_nr_from_lods(&marker[j],i,x2);
  413.           SHFT(f0,f1,f2,f)
  414.           } 
  415.         else
  416.           {
  417.           SHFT(x3,x2,x1,GRATIO*x2+MGRATIO*x0)
  418.           for (f=0.0,j=0;j<num_markers;++j)
  419.              if (marker[j].npairs[i]>2) f+=get_nr_from_lods(&marker[j],i,x1);
  420.           SHFT(f3,f2,f1,f)
  421.           }
  422.         }
  423.       if (f1 < f2)
  424.         { fmin=f1; minrel=x1; }
  425.       else
  426.         { fmin=f2; minrel=x2; }
  427. /* if best fitting value is close to 0.51 or to 0.99, check to see if 
  428. one of these values actually gives a better fit */
  429.       if (minrel<0.6)
  430.         {
  431.         for (f=0.0,j=0;j<num_markers;++j)
  432.              if (marker[j].npairs[i]>2) f+=get_nr_from_lods(&marker[j],i,0.51);
  433.         if (f<fmin) minrel=0.51;
  434.         }
  435.       else if (minrel>0.9)
  436.         {
  437.         for (f=0.0,j=0;j<num_markers;++j)
  438.              if (marker[j].npairs[i]>2) f+=get_nr_from_lods(&marker[j],i,0.99);
  439.         if (f<fmin) minrel=0.99;
  440.         }
  441.       if (minrel<0.51) minrel=0.51;
  442.       else if (minrel>0.99) minrel=0.99;
  443.       }
  444.     reliability[i]=minrel;
  445.     printf("Fitted \"reliability\" for pedigree %d: %5.3f\n",i+1,reliability[i]);
  446.     for (j=0;j<num_markers;++j)
  447.        get_nr_from_lods(&marker[j],i,reliability[i]);
  448.     }
  449.   if (debug_file)
  450.      {
  451.      fprintf(debug_file,"ped %3d, \"reliability\" = %5.3f:\n",
  452.        i+1,reliability[i]);
  453.      for (j=0;j<num_markers;++j)
  454.         fprintf(debug_file,"%s: %6.3f nonrec, %6.3f rec\n",
  455.           marker[j].name,marker[j].nonrec[i],marker[j].rec[i]);
  456.      }
  457.   }
  458. return 1;
  459. }
  460.  
  461. /* following routine needs to find the best fitting value for the 
  462. proportion of recombinants given the supplied lod scores, etc. */
  463.  
  464. float get_rf_fraction(pairs,npairs,rel,ssq,tol)
  465. struct lod_rf_t pairs[];
  466. int npairs;
  467. double rel,tol;
  468. float *ssq;
  469. {
  470. /* BEWARE - this curve has local minima!! */
  471. /* hence the following pathetic search */
  472. float f0,f1,f2,f3,x0,x1,x2,x3,fmin,xmin;
  473. int i;
  474. for (fmin=BIG,i=0;i<=10;++i)
  475.   {
  476. /* begin with finding the approximate position of the global minimum */
  477.   f0=err_rf_fraction(pairs,npairs,rel,i/10.0);
  478.   if (f0<fmin) { fmin=f0; xmin=i/10.0; }
  479.   }
  480. /* now set up starting points for golden ratio minimisation and use 
  481. that to find a close approximation */
  482. if (xmin==0/10.0)
  483.   { x0=0.0; x3=0.1; }
  484. else if (xmin==10/10.0) 
  485. /* Yes, I know this _should_ be 1, but let's play safe... */
  486.   { x0=0.9; x3=1.0; }
  487. else
  488.   { x0=xmin-0.1; x3=xmin+0.1; }
  489. x1=x0+(x3-x0)*MGRATIO;
  490. x2=x0+(x3-x0)*GRATIO;
  491. tol/=2.0; /* because I will use it as the distance between x1 and x2 */
  492. f1=err_rf_fraction(pairs,npairs,rel,x1);
  493. f2=err_rf_fraction(pairs,npairs,rel,x2);
  494. while (tol<fabs(x1-x2))
  495.   {
  496.   if (f2 < f1)
  497.     {
  498.     SHFT(x0,x1,x2,GRATIO*x1+MGRATIO*x3)
  499.     SHFT(f0,f1,f2,err_rf_fraction(pairs,npairs,rel,x2))
  500.     } 
  501.   else
  502.     {
  503.     SHFT(x3,x2,x1,GRATIO*x2+MGRATIO*x0)
  504.     SHFT(f3,f2,f1,err_rf_fraction(pairs,npairs,rel,x1))
  505.     }
  506.   }
  507. if (f1 < f2)
  508.   { fmin=f1; xmin=x1; }
  509. else
  510.   { fmin=f2; xmin=x2; }
  511. if (xmin<0.1  && (f0=err_rf_fraction(pairs,npairs,rel,0.0))<fmin)
  512.   { fmin=f0; xmin=0.0; }
  513. else if (xmin>0.9  && (f3=err_rf_fraction(pairs,npairs,rel,1.0))<fmin)
  514.   { fmin=f3; xmin=1.0; }
  515. *ssq=fmin;
  516. return xmin;
  517. }
  518.  
  519. /* calculate how closely the supplied proportion of recombinants (ratio) 
  520. fits to the observed lod scores, and return the error - the sum of 
  521. squares of distances between the observed lod scores and those produced
  522. using ratio (elod)*/
  523.  
  524. float err_rf_fraction(pairs,npairs,rel,ratio)
  525. struct lod_rf_t pairs[];
  526. int npairs;
  527. double rel,ratio;
  528. {
  529. float tot_lod,tot_elod,elod[MAXPAIRS],SSq;
  530. int i;
  531. SSq=0.0;
  532. for (tot_elod=tot_lod=0.0,i=0;i<npairs;++i)
  533.   {
  534.   tot_lod+=pairs[i].lod;
  535.   elod[i]=log_bin_like(ratio,1.0,pairs[i].rf*rel+(1-pairs[i].rf)*(1-rel));
  536.   elod[i]-=log_bin_like(ratio,1.0,0.5);
  537.   
  538.   tot_elod+=elod[i];
  539.   }
  540. tot_lod=fabs(tot_lod);
  541. tot_elod=fabs(tot_elod);
  542. for (i=0;i<npairs;++i)
  543.   {
  544.   float err;
  545.   err=pairs[i].lod/tot_lod-elod[i]/tot_elod;
  546.   SSq+=err*err;
  547.   }
  548. return SSq;
  549. }
  550.  
  551. /* given the observed lod scores and proportion of recombinants (ratio), 
  552. calculate the actual total number of informative meioses */
  553.  
  554. float get_num_meioses(pairs,npairs,rel,ratio)
  555. struct lod_rf_t pairs[];
  556. int npairs;
  557. double rel,ratio;
  558. {
  559. float sigma_eo,sigma_e2,elod;
  560. int i;
  561. for (sigma_eo=sigma_e2=0.0,i=0;i<npairs;++i)
  562.   {
  563.   elod=log_bin_like(ratio,1.0,pairs[i].rf*rel+(1-pairs[i].rf)*(1-rel));
  564.   elod-=log_bin_like(ratio,1.0,0.5);
  565.   sigma_e2+=elod*elod;        
  566.   sigma_eo+=elod*pairs[i].lod;
  567.   }
  568. return max(sigma_eo/sigma_e2,0.0);
  569. /* this is an analytic solution which minimises the least squares 
  570. distance between observed lod scores and those generated from supplied 
  571. values of rel and ratio */
  572. /* don't allow a negative number! */
  573. }
  574.  
  575. input_data()
  576. {
  577. int i,num_markers;
  578. if (!open_files()) return 0;
  579. if (!get_num_fams()) return 0;
  580. if (!input_disease_info(&disease)) return 0;
  581. if ((marker=(struct marker_t *)calloc(MAXMARKERS,sizeof(struct marker_t)))==NULL)
  582.       return error("Not enough memory for marker array","");
  583. for (i=0,num_markers=0;input_marker(&marker[i],num_peds);++i)
  584.         {
  585.         if (i && marker[i].position<marker[i-1].position)
  586.           return error("Marker is out of sequence: ",marker[i].name);
  587.         ++num_markers;
  588.         }
  589. return num_markers;
  590. }
  591.  
  592. open_files()
  593. {
  594. int nargs;
  595. char line[241],output_file_name[121],graph_file_name[121],debug_file_name[121];
  596. prompt("Enter output file name +/- graph file name\n");
  597. get_line(line,240);
  598. output_file_name[0]=graph_file_name[0]=debug_file_name[0]='\0';
  599. output_file=graph_file=debug_file=NULL;
  600. nargs=sscanf(line,"%120s %120s %120s",output_file_name,graph_file_name,debug_file_name);
  601. if (nargs<1)
  602.    return error("No file name given ",line);
  603. else if ((output_file=fopen(output_file_name,"w"))==0)
  604.    return error("Error trying to open file ",output_file_name);
  605. else if (nargs>1 && (graph_file=fopen(graph_file_name,"w"))==0)
  606.    return error("Error trying to open file ",graph_file_name);
  607. else if (nargs>2 && (debug_file=fopen(debug_file_name,"w"))==0)
  608.    return error("Error trying to open file ",debug_file_name);
  609. else return 1;
  610. }
  611.  
  612. set_up_tables(num_peds,num_points)
  613. int num_peds,num_points;
  614. {
  615. int i;
  616. if ((totlods=(float *)calloc(num_points,sizeof(float)))==NULL)
  617.       return error("Not enough memory for totlods","");
  618. else if ((tablods=(float **)calloc(num_peds,sizeof(float*)))==NULL)
  619.       return error("Not enough memory for tablods","");
  620. else for (i=0;i<num_peds;++i)
  621.       if ((tablods[i]=(float *)calloc(num_points,sizeof(float)))==NULL)
  622.           return error("Not enough memory for tablods","");
  623. return 1;
  624. }
  625.  
  626. input_disease_info(d)
  627. struct marker_t *d;
  628. {
  629. char line[501],buff[501],rest[501];
  630. float temp;
  631. int nargs,i;
  632. prompt(
  633. "Enter disease locus name and either one \"reliability\" value for all\n\
  634. pedigrees, or a \"reliability\" value for each pedigree, or no value:\n");
  635. get_line(line,240);
  636. nargs=sscanf(line,"%20s %f %[^\n]",d->name,&temp,rest);
  637. i=0;
  638. if (nargs<1)
  639.    return error("Bad argument format:\n",line);
  640. else if (nargs==1)
  641.    for (i=0;i<num_peds;++i) reliability[i]=-1.0;
  642.       /* calculate best-fitting value later */
  643. else 
  644.    {
  645.    do {
  646.    if (temp<0.0 || temp>1.0)
  647.       return error("\"Reliability\" must lie between 0 and 1:\n",line);
  648.    reliability[i++]=temp;
  649.    strcpy(buff,rest);
  650.    *rest='\0';
  651.    nargs=sscanf(buff," %f %[^\n]",&temp,rest);
  652.       } while (nargs>=1);
  653.    if (i==1)
  654.       for (i=1;i<num_peds;++i) reliability[i]=reliability[0];
  655.    else if (i!=num_peds)
  656.        return error("Wrong number of reliability values:\n",line);
  657.    }
  658. return 1;
  659. }
  660.  
  661. get_num_fams()
  662. {
  663. char line[501],buff[501],rest[501];
  664. int nargs;
  665. prompt("Enter minimum and maximum position:\n");
  666. get_line(line,500);
  667. if ((nargs=sscanf(line,"%f %f",&mind,&maxd))!=2)
  668.    return error("Bad initial data:\n",line);
  669. prompt("Enter number of points to use or values for fixed positions:\n");
  670. get_line(line,500);
  671. if ((nargs=sscanf(line,"%f %[^\n]",&fixdist[0],rest))<1)
  672.    return error("Bad initial data:\n",line);
  673. nfixdist=1;
  674. if (nargs>1) do {
  675. strcpy(buff,rest);
  676. if ((nargs=sscanf(buff,"%f %[^\n]",&fixdist[nfixdist],rest))<1) break;
  677. if (fixdist[nfixdist]<=fixdist[nfixdist-1])
  678.   return error("Fixed distances must be entered in ascending order:\n",line);
  679. else ++nfixdist;
  680. } while (nargs==2);
  681. if (nfixdist==1)
  682.   {
  683.   num_points=fixdist[0];
  684.   nfixdist=0;
  685.   }
  686. else num_points=nfixdist;
  687. prompt("Enter number of pedigrees (or 1 for just total lods):\n");
  688. get_line(line,240);
  689. num_peds=0;
  690. if (sscanf(line,"%d",&num_peds)!=1)
  691.    return error("No number of pedigrees given:\n",line);
  692. else if (num_peds<1)
  693.    return error("Number of pedigrees must be positive:\n",line);
  694. else return num_peds;
  695. }
  696.  
  697. /* currently all errors are fatal, maybe they could be trapped and 
  698. dealt with if desired */
  699.  
  700. error_func(l,f,s1,s2)
  701. unsigned l;
  702. char *f,*s1,*s2;
  703. {
  704. fprintf(stderr,"Error on line %d of source file %s:\n %s%s\n\n",l,f,s1,s2);
  705. exit(1);
  706. return 0;
  707. }
  708.  
  709. calculate_lods(mind,maxd,num_points,disease,marker,num_markers)
  710. double mind,maxd;
  711. int num_points;
  712. struct marker_t *disease,marker[];
  713. int num_markers;
  714. {
  715. int i;
  716. for (i=0;i<num_peds;++i)
  717.    {
  718.    total_meioses[i]=estimate_total_meioses(marker,num_markers,i);
  719.    make_map(mind,maxd,num_points,disease,marker,num_markers,i);
  720.    }
  721. return i;
  722. }
  723.  
  724. /* get the number of meioses informative for the disease locus in a
  725. pedigree, based on the number of meioses also informative for each
  726. marker and the probability each marker has of being informative
  727. */
  728.  
  729. float estimate_total_meioses(marker,num_markers,pednum)
  730. struct marker_t marker[];
  731. int num_markers,pednum;
  732. {
  733. float N,xtot,nmeioses;
  734. float a,t,sigma_a2,sigma_t2;
  735. int i;
  736. for (i=0;i<num_markers;++i)
  737.    if (marker[i].prob_inf>0.999) marker[i].prob_inf=0.999;
  738. /* minimise sigma{log(sq(Na-t)/Nab)} */
  739. /* i.e. based on normal approximation to binomial function */
  740. /* the solution comes out the same if you use a chi-squared function 
  741. to approximate the likelihood, i.e. if you seek to minimise:
  742. sigma{sqr(t-Na)/Na+sqr(N-t-Nb)/Nb} */
  743. for (sigma_a2=sigma_t2=0.0,i=0;i<num_markers;++i)
  744.    {
  745.    if (marker[i].npairs[pednum]==0) continue;
  746.    a=marker[i].prob_inf;
  747.    t=marker[i].rec[pednum]+marker[i].nonrec[pednum];
  748.    sigma_t2+=t*t/(a*(1-a));
  749.    sigma_a2+=a/(1-a);
  750. /* i.e.   sigma_a2+=a*a/(a*(1-a)); */
  751.    }
  752. nmeioses=sqrt(sigma_t2/sigma_a2);
  753. /* this is an analytic solution for the above minimisation */
  754. for (i=0;i<num_markers;++i)
  755.    nmeioses=max(nmeioses,(marker[i].rec[pednum]+marker[i].nonrec[pednum])/0.99);
  756. /* make sure that no pedigree has less informative meioses than any 
  757. total associated with any marker */
  758. if (debug_file)
  759.   fprintf(debug_file,"\nEstimated total informative meioses for ped %3d: %f\n",
  760.       pednum+1,nmeioses);
  761. for (i=0;i<num_markers;++i)
  762.    {
  763.    marker[i].fraction_used=
  764.      nmeioses ? (marker[i].rec[pednum]+marker[i].nonrec[pednum])/nmeioses : 0.0;
  765.    if (debug_file)
  766.        fprintf(debug_file,"%s.fraction_used: %f\n",
  767.            marker[i].name,marker[i].fraction_used);
  768.    }
  769. return nmeioses;
  770. }
  771.  
  772. /* initially fill and info2p structure with information from a marker 
  773. */
  774.  
  775. void info_from_marker(d,s,pednum)
  776. info2p_t *d;
  777. struct marker_t *s;
  778. int pednum;
  779. {
  780. d->position=s->position;
  781. d->fraction_used=s->fraction_used;
  782. d->rec=s->rec[pednum];
  783. d->nonrec=s->nonrec[pednum];
  784. }
  785.  
  786. /* here's a tricky bit - try to find a fair way of distributing the 
  787. meioses relating to two markers into nine possible categories 
  788. according to whether they are nonrecombinant, recombinant or 
  789. uninformative for each */
  790.  
  791. void sort_two_infos(l,r,result)
  792. info2p_t *l,*r;    
  793. float result[9];
  794. {
  795. float theta,tot_m;
  796. int i,j;
  797. theta=inf_rf_between(l,r);
  798. if (r->fraction_used==0.0 || l->fraction_used==0.0 || 
  799.     r->fraction_used*l->fraction_used<MINFRACTION) 
  800. {
  801. /* if only a small proportion of meioses are expected to be informative
  802. for both markers, ignore the overlap and pretend it is none at all */
  803. result[rr]=result[rn]=result[nr]=result[nn]=0.0;
  804. result[ru]=l->rec;
  805. result[nu]=l->nonrec;
  806. result[ur]=r->rec;
  807. result[un]=r->nonrec;
  808. if (r->fraction_used==0.0 && l->fraction_used==0.0)
  809.   result[uu]=0.0;
  810. else
  811.   {
  812.   result[uu]=r->fraction_used ? (r->rec+r->nonrec)/r->fraction_used
  813.                               : (l->rec+l->nonrec)/l->fraction_used;
  814.   for (i=0;i<8;++i) result[uu]-=result[i];
  815.   }
  816. }
  817. else
  818.   {
  819. /* select initial estimates for numbers falling into each category */
  820.   float expect[3][3],tot_r[3],tot_c[3];
  821.   int iter;
  822.   expect[0][0]= 2.0 / (1/(l->nonrec*(1-theta)*r->fraction_used) +
  823.                 1/(r->nonrec*(1-theta)*l->fraction_used));
  824.   expect[0][1]= 2.0 / (1/(l->nonrec*theta*r->fraction_used) +
  825.                 1/(r->rec*theta*l->fraction_used));
  826.   expect[1][0]= 2.0 / (1/(l->rec*theta*r->fraction_used) +
  827.                 1/(r->nonrec*theta*l->fraction_used));
  828.   expect[1][1]= 2.0 / (1/(l->rec*(1-theta)*r->fraction_used) +
  829.                 1/(r->rec*(1-theta)*l->fraction_used));
  830.   expect[0][2]=l->nonrec*(1-r->fraction_used);
  831.   expect[1][2]=l->rec*(1-r->fraction_used);
  832.   expect[2][0]=r->nonrec*(1-l->fraction_used);
  833.   expect[2][1]=r->rec*(1-l->fraction_used);
  834.   tot_m=(l->nonrec+l->rec)/l->fraction_used;
  835.   expect[2][2]=(1-l->fraction_used)*(1-r->fraction_used)*tot_m;
  836.  
  837. /* set up constraints row and column totals should conform to */
  838.   tot_r[0]=l->nonrec;
  839.   tot_r[1]=l->rec;
  840.   tot_r[2]=tot_m-l->nonrec-l->rec;
  841.   tot_c[0]=r->nonrec;
  842.   tot_c[1]=r->rec;
  843.   tot_c[2]=tot_m-r->nonrec-r->rec;
  844.   
  845. /* iterate starting values to fit to these constraints, making 
  846. proportional adjustments of the values on the worst-fitting row and 
  847. column */
  848. /* each column and row should contain at least one value to have got to 
  849. thsi stage */
  850. #define MAXERRTOL 0.01
  851. #define NITER 100
  852. /* seems to converge OK with these parameters */
  853.   for (iter=0;iter<NITER;++iter)
  854.     {
  855.     float total,maxerr;
  856.  
  857.     maxerr=-1;
  858.     for (total=0.0,i=0;i<3;++i)
  859.       {
  860.       if (tot_r[i]==0.0) continue;
  861.       total=expect[i][0]+expect[i][1]+expect[i][2];
  862.       if (total==0) error("Trying to adjust matrix with zero row total","");
  863.       maxerr=max(fabs(tot_r[i]-total)/tot_r[i],maxerr);
  864.       for (j=0;j<3;++j)
  865.         expect[i][j]*=tot_r[i]/total;
  866.       }
  867.     if (maxerr<MAXERRTOL) break;
  868.       maxerr=-1;
  869.       for (total=0.0,i=0;i<3;++i)
  870.         {
  871.         if (tot_c[i]==0.0) continue;
  872.         total=expect[0][i]+expect[1][i]+expect[2][i];
  873.         if (total==0) error("Trying to adjust matrix with zero column total","");
  874.         maxerr=max(fabs(tot_c[i]-total)/tot_c[i],maxerr);
  875.         for (j=0;j<3;++j)
  876.           expect[j][i]*=tot_c[i]/total;
  877.         }
  878.     if (maxerr<MAXERRTOL) break;
  879.     }
  880.   if (iter==NITER) error("Matrix adjustment failed to converge","");
  881.   result[nn]=expect[0][0];
  882.   result[nr]=expect[0][1];
  883.   result[nu]=expect[0][2];
  884.   result[rn]=expect[1][0];
  885.   result[rr]=expect[1][1];
  886.   result[ru]=expect[1][2];
  887.   result[un]=expect[2][0];
  888.   result[ur]=expect[2][1];
  889.   result[uu]=expect[2][2];
  890.   }
  891. }
  892.  
  893. /* adjust info structure so that the fraction of total meioses for 
  894. which this marker is deemed to be informative is correct again (since 
  895. rec and nonrec will have been adjusted to ignore information for those 
  896. meioses for which a closer marker was informative) */
  897.  
  898. void fix_fraction_used(i,tm)
  899. info2p_t *i;
  900. double tm;
  901. {
  902. i->fraction_used=min((i->rec+i->nonrec)/tm,1.0);
  903. }
  904.  
  905. void fill_info3p(p,l,r,rrnn)
  906. info3p_t *p;
  907. info2p_t *l,*r;
  908. float *rrnn;
  909. {
  910. int i;
  911. p->l_position=l->position;
  912. p->r_position=r->position;
  913. for (i=0;i<4;++i)
  914.    p->rrnn[i]=rrnn[i];
  915. }
  916.  
  917. /* go through calculating a lod score at each desired distance */
  918.  
  919. make_map(mind,maxd,num_points,disease,marker,num_markers,pednum)
  920. double mind,maxd;
  921. int num_points;
  922. struct marker_t *disease,marker[];
  923. int num_markers,pednum;
  924. {
  925. int i,j,k,new_interval,num_left,num_right;
  926. info2p_t l_info[MAXMARKERS],r_info[MAXMARKERS],tot_info;
  927. info3p_t info_table[MAXMARKERS][MAXMARKERS];
  928. new_interval=1;
  929. num_left=0;
  930. num_right=num_markers;
  931.  
  932. for (i=0,disease->position=nfixdist?fixdist[0]:mind;
  933.    i<num_points;
  934.    disease->position=nfixdist?fixdist[i+1]:
  935.    disease->position+(maxd-mind)/num_points,++i)
  936.   {
  937.   float lod;
  938. /* if disease position has moved past a marker then rearrange markers on
  939. either side and go on to sort out which meioses are informative for 
  940. which markers */
  941.   while (disease->position>marker[num_left].position && num_right)
  942.      {
  943.      ++num_left;
  944.      --num_right;
  945.      new_interval=1;
  946.      }
  947.   if (new_interval)
  948.     {
  949.     float result[9],theta;
  950.     new_interval=0;
  951.     for (j=0;j<num_left;++j)
  952.        info_from_marker(&l_info[j],&marker[num_left-j-1],pednum);
  953.     for (j=0;j<num_markers-num_left;++j)
  954.        info_from_marker(&r_info[j],&marker[num_left+j],pednum);
  955. /* begin by working through groups of markers to left and right of 
  956. disease locus separately */
  957.     if (num_left>0)
  958.       {
  959. /* tot_info relates to a "virtual marker" at the position of the last 
  960. marker and refers to the total number of meioses for which any closer 
  961. marker has been informative, and the expected state of these meioses (
  962. nonrecombinant or recombinant with the disease locus) */
  963.       tot_info.rec=l_info[0].rec;
  964.       tot_info.nonrec=l_info[0].nonrec;
  965.       tot_info.position=l_info[0].position;
  966.       tot_info.fraction_used=l_info[0].fraction_used;
  967.       for (j=1;j<num_left;++j)
  968.         {
  969.         sort_two_infos(&tot_info,&l_info[j],result);
  970.         theta=inf_rf_between(&tot_info,&l_info[j]);
  971. /* ignore meioses for which a closer marker has been informative */
  972.         l_info[j].rec=result[ur];
  973.         l_info[j].nonrec=result[un];
  974.         fix_fraction_used(&l_info[j],total_meioses[pednum]);
  975.         if (l_info[j].fraction_used<MINFRACTION)
  976.             clear_2p(&l_info[j]);
  977. /* if marker is only informative for a small proportion of meioses, pretend
  978. it is none */
  979. /* now update tot_info */
  980.         tot_info.rec=result[rr]+result[nr]+l_info[j].rec+
  981.            result[ru]*(1-theta)+result[nu]*theta;
  982.         tot_info.nonrec=result[nn]+result[rn]+l_info[j].nonrec+
  983.            result[nu]*(1-theta)+result[ru]*theta;
  984.         fix_fraction_used(&tot_info,total_meioses[pednum]);
  985.         }
  986.       }
  987.     if (num_right>0)
  988.       {
  989.       tot_info.rec=r_info[0].rec;
  990.       tot_info.nonrec=r_info[0].nonrec;
  991.       tot_info.position=r_info[0].position;
  992.       tot_info.fraction_used=r_info[0].fraction_used;
  993.       for (j=1;j<num_right;++j)
  994.         {
  995.         sort_two_infos(&tot_info,&r_info[j],result);
  996.         theta=inf_rf_between(&tot_info,&r_info[j]);
  997.         r_info[j].rec=result[ur];
  998.         r_info[j].nonrec=result[un];
  999.         fix_fraction_used(&r_info[j],total_meioses[pednum]);
  1000.         if (r_info[j].fraction_used<MINFRACTION)
  1001.             clear_2p(&r_info[j]);
  1002.             /* ignore this contribution */
  1003.         tot_info.rec=result[rr]+result[nr]+r_info[j].rec+
  1004.            result[ru]*(1-theta)+result[nu]*theta;
  1005.         tot_info.nonrec=result[nn]+result[rn]+r_info[j].nonrec+
  1006.            result[nu]*(1-theta)+result[ru]*theta;
  1007.         fix_fraction_used(&tot_info,total_meioses[pednum]);
  1008.         }
  1009.       }
  1010. /* now work through left and right groups together to find meioses for 
  1011. which markers on both sides of disease locus are informative */
  1012.     for (k=0;k<num_left;++k)
  1013.        {
  1014.        for (j=0;j<num_right;++j)
  1015.           {
  1016.           sort_two_infos(&l_info[k],&r_info[j],result);
  1017.           fill_info3p(&info_table[k][j],&l_info[k],&r_info[j],result);
  1018. /* discard meioses informative for both markers from further 
  1019. consideration before going on to more distant markers */
  1020.           l_info[k].rec=result[ru];
  1021.           l_info[k].nonrec=result[nu];
  1022.           fix_fraction_used(&l_info[k],total_meioses[pednum]);
  1023.           r_info[j].rec=result[ur];
  1024.           r_info[j].nonrec=result[un];
  1025.           fix_fraction_used(&r_info[j],total_meioses[pednum]);
  1026.           }
  1027.        }
  1028. if (debug_file)
  1029. {
  1030. fprintf(debug_file,"\n\n\n%8.8s "," ");
  1031. for(j=0;j<num_right;++j)
  1032.     fprintf(debug_file,"%8.3f  ",r_info[j].nonrec);
  1033. fprintf(debug_file,"\n%8.8s "," ");
  1034. for(j=0;j<num_right;++j)
  1035.     fprintf(debug_file,"  %8.3f",r_info[j].rec);
  1036. for (k=0;k<num_left;++k)
  1037.   {
  1038.   fprintf(debug_file,"\n\n%8.3f ",l_info[k].nonrec);
  1039.   for(j=0;j<num_right;++j)
  1040.       fprintf(debug_file,"%8.3f  ",info_table[k][j].rrnn[nn]);
  1041.   fprintf(debug_file,"\n%8.8s "," ");
  1042.   for(j=0;j<num_right;++j)
  1043.       fprintf(debug_file,"  %8.3f",info_table[k][j].rrnn[nr]);
  1044.   fprintf(debug_file,"\n%8.3f ",l_info[k].rec);
  1045.   for(j=0;j<num_right;++j)
  1046.       fprintf(debug_file,"%8.3f  ",info_table[k][j].rrnn[rn]);
  1047.   fprintf(debug_file,"\n%8.8s "," ");
  1048.   for(j=0;j<num_right;++j)
  1049.       fprintf(debug_file,"  %8.3f",info_table[k][j].rrnn[rr]);
  1050.   }
  1051. }
  1052. /* all meioses categorised now, go on to calculate lod scores */
  1053.   for (k=0;k<num_left;++k)
  1054.     set_null_like_2p(&l_info[k]);
  1055.   for(j=0;j<num_right;++j)
  1056.     set_null_like_2p(&r_info[j]);
  1057.   for (k=0;k<num_left;++k)
  1058.      for(j=0;j<num_right;++j)
  1059.        set_null_like_3p(&info_table[k][j]);
  1060.     }
  1061.   lod=0.0;
  1062.   for (k=0;k<num_left;++k)
  1063.     lod+=lod_2p(disease,&l_info[k],reliability[pednum]);
  1064.   for(j=0;j<num_right;++j)
  1065.     lod+=lod_2p(disease,&r_info[j],reliability[pednum]);
  1066.   for (k=0;k<num_left;++k)
  1067.      for(j=0;j<num_right;++j)
  1068.        lod+=lod_3p(disease,&info_table[k][j],reliability[pednum]);
  1069.   tablods[pednum][i]=max(lod,-99.0);
  1070.   }
  1071. return 1;
  1072. }
  1073.  
  1074. float log_bin_like(x,N,a)
  1075. double x,N,a;
  1076. {
  1077. return x*(a>0.0?log10(a):-BIG)+(N-x)*(1-a>0.0?log10(1-a):-BIG);
  1078. }
  1079.  
  1080. void prompt(s)
  1081. char *s;
  1082. {
  1083. if (input_file==stdin)
  1084.    printf("%s",s);
  1085. }
  1086.  
  1087. /* write table and optionally graph file */
  1088. /* this function could easily be adapted to produce output in whatever 
  1089. form was convenient */
  1090.  
  1091. void output_tables(num_peds,num_points,num_markers,mind,maxd)
  1092. int num_peds,num_points,num_markers;
  1093. double mind,maxd;
  1094. {
  1095. int ped,point,i;
  1096. float min_lod,max_lod;
  1097. for (point=0;point<num_points;++point)
  1098.   {
  1099.   totlods[point]=0.0; /* should be already */
  1100.   for (ped=0;ped<num_peds;++ped)
  1101.     totlods[point]+=tablods[ped][point];
  1102.   }
  1103. fprintf(output_file,"%-9s  %-8s","Distance",num_peds>1?"Total":"Lod");
  1104. if (num_peds>1)
  1105.    for (ped=0;ped<num_peds;++ped) fprintf(output_file," Ped %3d",ped+1);
  1106. for (point=0;point<num_points;++point)
  1107.   {
  1108.   fprintf(output_file,"\n%9.3f  %8.4f",
  1109.      nfixdist?fixdist[point]:mind+point*(maxd-mind)/num_points,
  1110.      totlods[point]);
  1111.   if (num_peds>1)
  1112.     for (ped=0;ped<num_peds;++ped) fprintf(output_file," %7.4f",tablods[ped][point]);
  1113.   }
  1114. fclose(output_file);
  1115. if (!graph_file) return;
  1116. fprintf(graph_file,"PARMS:L\n\
  1117. TITLE:FASTMAP OUTPUT: LINKAGE MAP FOR %s\n\
  1118. TITLEV:lod score\nTITLEG:cM\nTITLEG:cM\nTITLEG:%s",
  1119.   disease.name,
  1120.   num_peds>1?"Total":"Lod");
  1121. if (num_peds>1) for (ped=0;ped<num_peds;++ped)
  1122.    fprintf(graph_file,"\nTITLEG:ped %d",ped+1);
  1123. min_lod=max_lod=0.0;
  1124. for (point=0;point<num_points;++point)
  1125.   {
  1126.   fprintf(graph_file,"\n%f,%f",
  1127.     nfixdist?fixdist[point]:mind+point*(maxd-mind)/num_points,
  1128.     totlods[point]);
  1129.   min_lod=min(min_lod,totlods[point]);
  1130.   max_lod=max(max_lod,totlods[point]);
  1131.   if (num_peds>1)
  1132.     for (ped=0;ped<num_peds;++ped)
  1133.       {
  1134.       fprintf(graph_file,",%f",tablods[ped][point]);
  1135.       min_lod=min(min_lod,tablods[ped][point]);
  1136.       max_lod=max(max_lod,tablods[ped][point]);
  1137.       }
  1138.   }
  1139. fprintf(graph_file,"\nDATATYPE:XY1,2,-1");
  1140. if (num_peds>1)
  1141.   for (ped=0;ped<num_peds;++ped)
  1142.        fprintf(graph_file,"\nDATATYPE:XY1,%d,-%d",ped+3,ped+2);
  1143. fprintf(graph_file,"\nMINX:%.3f\nMAXX:%.3f",mind,maxd);
  1144. #if 0
  1145. if (min_lod<-3.0)
  1146.   {
  1147.   fprintf(graph_file,"\nMINY:-3.0");
  1148.   fprintf(graph_file,"\nMAXY:%.1f",ceil(max_lod));
  1149.   }
  1150. #endif
  1151. for (i=0;i<num_markers;++i)
  1152.   fprintf(graph_file,"\nTITLEF:%f,0.0,90,%s_",
  1153.      (marker[i].position-mind)/(maxd-mind),
  1154.      marker[i].name);
  1155. fclose(graph_file);
  1156. }
  1157.  
  1158. #ifndef TEST
  1159. /* allow testing of individual functions by other modules if desired */
  1160. main(argc,argv)
  1161. int argc;
  1162. char *argv[];
  1163. {
  1164. int num_markers;
  1165. print_consts();
  1166. if (argc>1)
  1167.    {
  1168.    input_file=fopen(argv[1],"r");
  1169.    if (!input_file)
  1170.      { error("Could not open file ",argv[1]); return 1; }
  1171.    }
  1172. else
  1173.   input_file=stdin;
  1174. if ((num_markers=input_data())!=0 &&
  1175. meioses_from_lods(num_peds,num_markers) &&
  1176. set_up_tables(num_peds,num_points) &&
  1177. calculate_lods(mind,maxd,num_points,&disease,marker,num_markers))
  1178.   output_tables(num_peds,num_points,num_markers,mind,maxd);
  1179. free(marker);
  1180. if (input_file!=stdin) fclose(input_file);
  1181. return 0;
  1182. }
  1183. #endif
  1184.  
  1185. 
  1186.